#How Rebels Keep the Peace
#Margaret McWeeney
#June 30, 2021

rm(list=objects()) 
library(MASS)
library(mediation)
library(Hmisc)
library(dplyr)
library(tidyverse)
library(stargazer)
library(psych)
library(effects)
library(survival)
library(ggplot2)

LegitConc<- read.csv("/Users/PeggyMcweeney/Desktop/LegitConcactiveV2_3_2_20.csv", header=TRUE)

LegitConcess2<-LegitConc

##Adding Time * time2 * Time3

LegitConcess2$time<-0
cut<-as.data.frame(NULL) #create empty data for each facid
res<-as.data.frame(NULL) #create empty data for the merged results
for(i in 1:length(unique(LegitConcess2$facid))){ #for each facid
  cut<-LegitConcess2[which(LegitConcess2$facid==unique(LegitConcess2$facid)[i]),] #subset the data for that facid only
  concessimp<-which(cut$concessimp==1) #which rows in this facid cut have terr?
  if(length(concessimp)>0){ #are there any terr events? if so,
    for (j in 1:length(concessimp)){ #for each event
      if(j==length(concessimp)){ #is this the last terrorist event for that facid? 
        #if so, from the year after that event, count to the end of the dataset
        #can't do it if the last year of the dataset, so only if it isn't, then...
        if(cut$year[concessimp[j]]!=max(cut$year)){
          cut$time[(concessimp[j]+1):dim(cut)[1]]<-1:length((concessimp[j]+1):dim(cut)[1])
        } #close if not the lat year of the dataset
      } #close if for last event
      #if it's not the last terrorist event, then 
      #from the year after that event through the year with the next one
      #insert numbers from 1 counting up to the number of years until the next event
      else{
        cut$time[(concessimp[j]+1):concessimp[j+1]]<-1:length((concessimp[j]+1):concessimp[j+1])
      } #close else
    } #close for each event
  } #close if for no terrorism
  res<-rbind(res, cut)  #bind this facid's results to the previous facids' results
} #close facid loop
LegitConcess2<-res
LegitConcess2$time<-LegitConcess2$time
LegitConcess2$time2<-LegitConcess2$time^2
LegitConcess2$time3<-LegitConcess2$time^3
#Use time, time^2, time^3 because concessions has an effect and then reduces over time if not used
#Restarts at the next incidence of concession so this has the smoothing effect
#For other variables I use regular lags to test whether they might have this effect

#Post-concession data
LegitConcess2<-subset(LegitConcess2, time>0)
LegitConcess2<-LegitConcess2[!(LegitConcess2$year=="1993"),]

#Descriptive Graphs
country<-as.data.frame(unique(LegitConcess2$country))

LegitConcess2$region<-NA
#Africa
LegitConcess2$region[LegitConcess2$country=="Chad"]<-"Africa"
LegitConcess2$region[LegitConcess2$country=="Comoros"]<-"Africa"
LegitConcess2$region[LegitConcess2$country=="Djibouti"]<-"Africa"
LegitConcess2$region[LegitConcess2$country=="Ethiopia"]<-"Africa"
LegitConcess2$region[LegitConcess2$country=="Somalia"]<-"Africa"
LegitConcess2$region[LegitConcess2$country=="Sudan"]<-"Africa"
LegitConcess2$region[LegitConcess2$country=="South Sudan"]<-"Africa"


#Central Asia
LegitConcess2$region[LegitConcess2$country=="Afghanistan"]<- "Central Asia"
LegitConcess2$region[LegitConcess2$country=="Pakistan"]<- "Central Asia"
LegitConcess2$region[LegitConcess2$country=="Russia"]<- "Central Asia"
LegitConcess2$region[LegitConcess2$country=="Georgia"]<- "Central Asia"

#East Asia

#Europe
LegitConcess2$region[LegitConcess2$country=="Macedonia"]<- "Europe"
LegitConcess2$region[LegitConcess2$country=="Spain"]<-"Europe"
LegitConcess2$region[LegitConcess2$country=="United Kingdom"]<-"Europe"
LegitConcess2$region[LegitConcess2$country=="Serbia"]<-"Europe"
LegitConcess2$region[LegitConcess2$country=="Moldova"]<-"Europe"

#Middle East
LegitConcess2$region[LegitConcess2$country=="Iraq"]<-"Middle East"
LegitConcess2$region[LegitConcess2$country=="Israel"]<-"Middle East"

#North America
LegitConcess2$region[LegitConcess2$country=="United States"]<-"North America"

#Oceania

#Southeast/South Asia
LegitConcess2$region[LegitConcess2$country=="Sri Lanka"]<-"South/S.east Asia"
LegitConcess2$region[LegitConcess2$country=="Bangladesh"]<-"South/S.east Asia"
LegitConcess2$region[LegitConcess2$country=="Myanmar"]<-"South/S.east Asia"
LegitConcess2$region[LegitConcess2$country=="Thailand"]<-"South/S.east Asia"
LegitConcess2$region[LegitConcess2$country=="Philippines"]<-"South/S.east Asia"
LegitConcess2$region[LegitConcess2$country=="Indonesia"]<-"South/S.east Asia"

#India
LegitConcess2$region[LegitConcess2$country=="India"]<-"India"

TRYING<-distinct(LegitConcess2,name, .keep_all= TRUE)
Orgcounts <- table(TRYING$region)

#Figure 1
Orgcounts<-as.data.frame(Orgcounts)
barplot(Orgcounts,
        main="Figure 1: Organizations by Region",
        xlab="World Regions",
        ylab="Number of Organizations",
        las=1,
        cex.names=0.50,
        names.arg=c("Africa", "Central Asia", "Europe", "India", "Middle East", "North America", "South/S.east Asia"))


p<-ggplot(Orgcounts, aes(x=Var1, y=Freq, fill=Var1)) + 
  geom_bar(stat = "identity")

p + xlab("World Region") + ylab("Number of Organizations") +
  guides(fill=FALSE, color=FALSE)
  

#Figure 2
Africa<-subset(TRYING, TRYING$region=="Africa")
mean(Africa$age)
CentralAsia<-subset(TRYING, TRYING$region=="Central Asia")
mean(CentralAsia$age)
Europe<-subset(TRYING, TRYING$region=="Europe")
mean(Europe$age)
MiddleEast<-subset(TRYING, TRYING$region=="Middle East")
mean(MiddleEast$age)
NorthAm<-subset(TRYING, TRYING$region=="North America")
mean(NorthAm$age)
SEAsia<-subset(TRYING, TRYING$region=="South/S.east Asia")
mean(SEAsia$age)
India<-subset(TRYING, TRYING$region=="India")
mean(India$age)

Region<-c("Africa", "Central Asia", "Europe", "India", "Middle East", "North America", "South/S.east Asia")
Meanage<-c(22.9, 24.7, 31, 38.2, 43.4, 15, 39.43)
TRYINGHARD <- data.frame(Region, Meanage)

barplot(TRYINGHARD$Meanage,
        main="Figure 2: Mean Organizational Age by Region",
        xlab="World Regions",
        ylab="Mean Age of Organizations",
        las=1,
        cex.names=0.50,
        names.arg=c("Africa", "Central Asia", "Europe", "India", "Middle East", "North America", "South/S.east Asia"))

#No plot - used for information in paper

#Post-Concession Terrorism

LegitConcess2$postconterrciv<-LegitConcess2$civtarg

LegitConcess2$postconterrstate<-LegitConcess2$statetarg

LegitConcess2$postcon_nkill<-LegitConcess2$nkill

LegitConcess2$postcon_nkillciv<-LegitConcess2$nkill
#LegitConcess2$postcon_nkillciv[LegitConcess2$postconterrciv>1]<-0

LegitConcess2$postcon_nkillstate<-LegitConcess2$nkill
#LegitConcess2$postcon_nkillstate[LegitConcess2$postconterrstate>1]<-0

## Variable lags

LegitConcess2<-LegitConcess2 %>% 
  group_by(name) %>% 
  mutate(terr5 = dplyr::lag(gtdterrevent, n=5, default=0))
LegitConcess2<-LegitConcess2 %>% 
  group_by(name) %>% 
  mutate(terr1 = dplyr::lag(gtdterrevent, n=1, default=0))
LegitConcess2<-LegitConcess2 %>% 
  group_by(name) %>% 
  mutate(terr3 = dplyr::lag(gtdterrevent, n=3, default=0))

#Figure 3
terrorism<-subset(LegitConcess2, LegitConcess2$gtdterrevent>0)
sum(terrorism$gtdterrevent)
myvars<-c("name", "country", "facid", "year", "region", "gtdterrevent")
terrorism<-terrorism[myvars]
terrcounts<-LegitConcess2 %>% group_by(region) %>% summarise(gtdterrevent = sum(gtdterrevent))
sum(terrcounts$gtdterrevent)
terrcounts <- terrcounts[-c(8), ]
terrcounts<-as.data.frame(terrcounts)
barplot(terrcounts$gtdterrevent,
        main="Figure 3: Post-concession Terrorism by Region",
        xlab="World Regions",
        ylab="Frequency of Post-concession Terrorism Events",
        las=1,
        cex.names=0.1,
        xaxs="i")

p<-ggplot(terrcounts, aes(x=region, y=gtdterrevent, fill=region)) + 
  geom_bar(stat = "identity")

p + xlab("World Region") + ylab("Frequency of Post-Concession Terrorism") +
  guides(fill=FALSE, color=FALSE)

length(unique(terrorism$name))
#42
##Middle East orgs
MEterrorism<-subset(terrorism, region=="Middle East")
length(unique(MEterrorism$name))
length(MEterrorism$gtdterrevent)

#South Asia
SEterrorism<-subset(terrorism, region=="South/S.east Asia")
length(unique(SEterrorism$name))
length(SEterrorism$gtdterrevent)

#Africa
Aterrorism<-subset(terrorism, region=="Africa")
length(unique(Aterrorism$name))
length(Aterrorism$gtdterrevent)

#Central Asia
CAterrorism<-subset(terrorism, region=="Central Asia")
length(unique(CAterrorism$name))
length(CAterrorism$gtdterrevent)

#Europe
Eterrorism<-subset(terrorism, region=="Europe")
length(unique(Eterrorism$name))
length(Eterrorism$gtdterrevent)

#India
Iterrorism<-subset(terrorism, region=="India")
length(unique(Iterrorism$name))
length(Iterrorism$gtdterrevent)

#Figure 4: Women in post-concession orgs
women<-subset(LegitConcess2, wom==1)
unique(women$name)
length(unique(women$name))

womencount<-distinct(women,name, .keep_all= TRUE)
womencount <- table(womencount$region)

womencount<-as.data.frame(womencount)

p<-ggplot(womencount, aes(x=Var1, y=Freq, fill=Var1)) + 
  geom_bar(stat = "identity")

p + xlab("World Region") + ylab("Count of Organizations Including Women") +
  guides(fill=FALSE, color=FALSE)

#Table 1, bivariate relationship between women and terrorism
length(unique(LegitConcess2$facid))
womenterr<-LegitConcess2 %>% group_by(facid) %>% summarise(gtdterrevent = sum(gtdterrevent))
womenterr1<-LegitConcess2 %>% group_by(facid) %>% summarise(nkill = sum(nkill))
womenterr2<-LegitConcess2 %>% group_by(facid) %>% summarise(wom = sum(wom))
womenterr3<-LegitConcess2 %>% group_by(facid) %>% summarise(rebelections = sum(rebelections))
womenterr4<-merge(womenterr, womenterr1, by="facid")
womenterr4<-merge(womenterr4, womenterr2, by="facid")
womenterr5<-merge(womenterr4, womenterr3, by="facid")


#Chi-square
womenterr5<- na.omit(womenterr5)
table(womenterr5$gtdterrevent, womenterr5$wom)
womenterr5$gtdterrevent[womenterr5$gtdterrevent>0] <- 1
tbl<-table(womenterr5$gtdterrevent, womenterr5$wom)
chisq.test(tbl)

#Figure 4: Rebel elections
LegitConcess2$rebelections2<-LegitConcess2$rebelections
LegitConcess2$rebelections2[LegitConcess2$rebelections>1]<-1
rebelections2<-subset(LegitConcess2, rebelections2==1)
unique(LegitConcess2$rebelections2)

rebelectionscount<-distinct(rebelections2,name, .keep_all= TRUE)
rebelectionscount <- table(rebelectionscount$region)

rebelectionscount<-as.data.frame(rebelectionscount)

p<-ggplot(rebelectionscount, aes(x=Var1, y=Freq, fill=Var1)) + 
  geom_bar(stat = "identity")

p + xlab("World Region") + ylab("Count of Organizations with Elections") +
  guides(fill=FALSE, color=FALSE)

##Appendix A - Descriptive Statistics

keeps <- c("autonomy", "armedconflict", "v2x_polyarchy", "gdp_pc",  "islamist",
           "gtdterrevent", "postconterrciv", "postconterrstate",
            "rebelections2", "wom")
LegitConcess3 <- LegitConcess2[keeps]
 
# write.csv(LegitConcess3,"/Users/PeggyMcweeney/Desktop/DSPaper3.csv")

###Count Model, clustered by organization
LegitConcess2$rebel_elections<-LegitConcess2$rebelections
LegitConcess2$rebel_elections[LegitConcess2$rebelections>1]<-1

install.packages('R2admb')
install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")
library(R2admb)
library(glmmADMB)
m1 <- glmmadmb(gtdterrevent ~ wom + rebel_elections
                     + wom:rebel_elections + islamist
                     + autonomy, 
                     data=LegitConcess2, zeroInflation=TRUE, family="poisson")
summary(m1)
exp(coef(m1))

m2 <- glmmadmb(gtdterrevent ~ wom + rebel_elections
               + wom:rebel_elections + islamist
               + autonomy
               + armedconflict + v2x_polyarchy 
               + gdp_pc  + time + time2 + time3, 
               data=LegitConcess2, zeroInflation=TRUE, family="poisson")

summary(m2)

exp(coef(m2))

##Appendix B - Fatalities instead of terrorism
m3 <- glmmadmb(nkill ~ wom + rebel_elections
               + wom:rebel_elections + islamist
               + autonomy, 
               data=LegitConcess2, zeroInflation=TRUE, family="poisson")
summary(m3)
exp(coef(m3))

m4 <- glmmadmb(nkill ~ wom + rebel_elections
               + wom:rebel_elections + islamist
               + autonomy
               + armedconflict + v2x_polyarchy 
               + gdp_pc  + time + time2 + time3, 
               data=LegitConcess2, zeroInflation=TRUE, family="poisson")

summary(m4)
exp(coef(m4))

#did not work - manual input into table
stargazer(m1,m2, title="Dependent variable:", align=TRUE, type="html", out="Paper3.html")

#Survival Model Restarting from beginning - re-entering base datafile
library("survival")
library("tidyverse")

LegitConc<- read.csv("/Users/PeggyMcweeney/Desktop/LegitConcactiveV2_3_2_20.csv", header=TRUE)

LegitConcess2<-LegitConc

##Adding Time * time2 * Time3

LegitConcess2$time<-0
cut<-as.data.frame(NULL) #create empty data for each facid
res<-as.data.frame(NULL) #create empty data for the merged results
for(i in 1:length(unique(LegitConcess2$facid))){ #for each facid
  cut<-LegitConcess2[which(LegitConcess2$facid==unique(LegitConcess2$facid)[i]),] #subset the data for that facid only
  concessimp<-which(cut$concessimp==1) #which rows in this facid cut have terr?
  if(length(concessimp)>0){ #are there any terr events? if so,
    for (j in 1:length(concessimp)){ #for each event
      if(j==length(concessimp)){ #is this the last terrorist event for that facid? 
        #if so, from the year after that event, count to the end of the dataset
        #can't do it if the last year of the dataset, so only if it isn't, then...
        if(cut$year[concessimp[j]]!=max(cut$year)){
          cut$time[(concessimp[j]+1):dim(cut)[1]]<-1:length((concessimp[j]+1):dim(cut)[1])
        } #close if not the lat year of the dataset
      } #close if for last event
      #if it's not the last terrorist event, then 
      #from the year after that event through the year with the next one
      #insert numbers from 1 counting up to the number of years until the next event
      else{
        cut$time[(concessimp[j]+1):concessimp[j+1]]<-1:length((concessimp[j]+1):concessimp[j+1])
      } #close else
    } #close for each event
  } #close if for no terrorism
  res<-rbind(res, cut)  #bind this facid's results to the previous facids' results
} #close facid loop
LegitConcess2<-res
LegitConcess2$time<-LegitConcess2$time
LegitConcess2$time2<-LegitConcess2$time^2
LegitConcess2$time3<-LegitConcess2$time^3
#Use time, time^2, time^3 because concessions has an effect and then reduces over time if not used
#Restarts at the next incidence of concession so this has the smoothing effect
#For other variables I use regular lags to test whether they might have this effect

#Post-concession data
LegitConcess2<-subset(LegitConcess2, time>0)
LegitConcess2<-LegitConcess2[!(LegitConcess2$year=="1993"),]

#Post-Concession Terrorism

LegitConcess2$postconterrciv<-LegitConcess2$civtarg

LegitConcess2$postconterrstate<-LegitConcess2$statetarg

LegitConcess2$postcon_nkill<-LegitConcess2$nkill

LegitConcess2$postcon_nkillciv<-LegitConcess2$nkill
#LegitConcess2$postcon_nkillciv[LegitConcess2$postconterrciv>1]<-0

LegitConcess2$postcon_nkillstate<-LegitConcess2$nkill
#LegitConcess2$postcon_nkillstate[LegitConcess2$postconterrstate>1]<-0

LegitConcess2<-LegitConcess2 %>% 
  group_by(name) %>% 
  mutate(terr5 = dplyr::lag(gtdterrevent, n=5, default=0))
LegitConcess2<-LegitConcess2 %>% 
  group_by(name) %>% 
  mutate(terr1 = dplyr::lag(gtdterrevent, n=1, default=0))
LegitConcess2<-LegitConcess2 %>% 
  group_by(name) %>% 
  mutate(terr3 = dplyr::lag(gtdterrevent, n=3, default=0))

#interaction variable
LegitConcess2$rebel_elections<-LegitConcess2$rebelections
LegitConcess2$rebel_elections[LegitConcess2$rebelections>1]<-1
LegitConcess2$interaction<-NA
LegitConcess2$interaction<-LegitConcess2$wom*LegitConcess2$rebel_elections

#######################Survival Analysis 
#https://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/R/R7_LogisticRegression-Survival/R7_LogisticRegression-Survival4.html

LegitConcess2 <- with(LegitConcess2, LegitConcess2[complete.cases(gtdterrevent, islamist, autonomy, terr1, terr3, terr5, 
                                                                  wom, rebel_elections, interaction, armedconflict, v2x_polyarchy, gdp_pc, time), ])
#postconterrstate
#postconterrciv
#islamist
#autonomy
#terr1
#terr3
#terr5
hist(LegitConcess2$time)
KM_fit <- survfit(Surv(LegitConcess2$time, LegitConcess2$gtdterrevent) ~ 1)
KM_fit
plot(KM_fit, xlab="Time to first concession (years)")
#Median=19 years to reach 50-50 chance of terrorism

#Baseline model, without interaction
coxph(formula=Surv(time, gtdterrevent>0)~ wom + rebel_elections
      + armedconflict + v2x_polyarchy + islamist
      + autonomy + terr1 + terr3 + terr5
      + gdp_pc + interaction)

m10 = coxph(Surv(LegitConcess2$time,LegitConcess2$gtdterrevent>0) ~ wom + rebel_elections
            + armedconflict + v2x_polyarchy + islamist
            + autonomy + terr1 + terr3 + terr5
            + gdp_pc, na.action = na.omit, data=LegitConcess2)
summary(m10)

Cox_curve1 <- survfit(m10)
plot(Cox_curve1, main="Figure 7: Proportional Hazard Rate of Terrorism Post-Concession",
     xlab="Years since concession", ylab="Probability of Terrorism",
     conf.int=FALSE, lty=1, xlim=c(0,10))
legend("bottomleft", legend=c("Baseline", "Women x Rebel Elections"),
       lty=c(1,3), box.lty=1) 

#Figure 7 - Plot with women x elections interaction
m11 = coxph(Surv(LegitConcess2$time,LegitConcess2$gtdterrevent>0) ~ wom + rebel_elections
            + armedconflict + v2x_polyarchy + islamist
            + autonomy + terr1 + terr3 + terr5 + interaction
            + gdp_pc, na.action = na.omit, data=LegitConcess2)
summary(m11)
Cox_curve2 <- survfit(m11)
par(new=TRUE)
plot(Cox_curve2, lty=3, frame.plot=FALSE, axes=FALSE, conf.int=FALSE, xlab="Time in years")

stargazer(m10,m11, title="Dependent variable:", align=TRUE, type="html", out="Paper3.1.html")

